#ifndef CALCULATION
#define CALCULATION
#include "basic.h"
#include <map>
#include <functional>
class data_calculation : public basic_valiable
{
public:
    double f_p_jins(double p);
    double px_py_p_chuiz(double p_zhuiz);
    double f_theta_1(double theta);
    double f_theta_2(double theta);
    double f_theta_3(double theta);
    double f_p(double p);
    double f_p_2(double p);
    double f_theta_integra(double begin, double end, double theta);
    double f_p_integra(double begin, double end, double p);
    void plot_px_py_p_chuiz();
};
double data_calculation::f_p_jins(double p)
{
    double sigma_theta = sigma_p_chuiz / p_0;
    double temp = pow(sigma_theta * p, 2) + pow(sigma_p_pingx, 2);
    double x_1 = exp(-(p_0 * p_0) / 2 / temp);
    double x_2 = pow(sigma_theta * p * p_0, 2) * sigma_theta + pow(sigma_theta * sigma_p_pingx * p, 2) * sigma_theta + pow(sigma_p_pingx, 4) * sigma_theta;
    return x_1 * x_2 / pow(temp, 2.5);
}
double data_calculation::f_theta_3(double theta)
{
    double sigma = sigma_p_chuiz / p_0;
    return theta / pow(sigma, 2) * exp(-pow(theta, 2) / 2 / pow(sigma, 2));
}
double data_calculation::f_theta_2(double theta)
{
    double x1 = tan(theta) * (pow(sigma_p_chuiz * p_0, 2) * sigma_p_chuiz + sigma_p_chuiz * pow(sigma_p_pingx, 4) * pow(tan(theta), 2) + pow(sigma_p_chuiz * sigma_p_pingx, 2) * sigma_p_chuiz);
    double temp = pow(sigma_p_chuiz, 2) + pow(sigma_p_pingx * tan(theta), 2);
    double x2 = pow(cos(theta), 2) * pow(temp, 2.5);
    double x3 = pow(p_0 * tan(theta), 2) / 2 / temp;
    return x1 / x2 * exp(-x3);
}
double data_calculation::f_p_2(double p)
{
    double up = 0.125;
    double sigma_theta = sigma_p_chuiz / p_0;
    double k = (pow(sigma_theta * p, 2) + pow(sigma_p_chuiz, 2)) / (2 * pow(sigma_p_chuiz * sigma_theta, 2));
    return p / 2 / pow(sigma_p_chuiz * sigma_theta, 2) * (1 / k / k - (up * up / k + 1 / k / k) * exp(-k * up * up));
}
double data_calculation::f_p_integra(double begin, double end, double p)
{
    double sum = 0;

    for (double i = begin; i < end; i += 1)
    {

        sum += exp(-((pow(p, 2) - pow(i, 2))) / (2 * pow(sigma_p_chuiz, 2))) * exp(-(pow(i - p_0, 2) / (2 * pow(sigma_p_pingx, 2))));
    }
    return sum;
}
double data_calculation::f_p(double p)
{
    return p / sqrt(2 * M_PI) / sigma_p_pingx / pow(sigma_p_chuiz, 2) * f_p_integra(-p, p, p);
}
double data_calculation::f_theta_integra(double begin, double end, double theta)
{
    double sum = 0;
    double temp = pow(sigma_p_pingx * tan(theta), 2) + pow(sigma_p_chuiz, 2);
    double mu = pow(sigma_p_chuiz, 2) * p_0 / temp;
    double sigma = pow(sigma_p_chuiz * sigma_p_pingx, 2) / temp;
    for (double p = begin; p < end; p += 0.1)
    {

        sum += 0.1 * pow(p, 2) * exp(-pow(p - mu, 2) / 2 / sigma);
    }
    return sum;
}
double data_calculation::f_theta_1(double theta)
{
    // return exp(-(pow(p_0 * tan(theta), 2) / 2 / (pow(sigma_p_pingx * tan(theta), 2) + pow(sigma_p_chuiz, 2)))) * sin(theta) / pow(cos(theta), 3) * (pow(sigma_p_chuiz, 3) * pow(p_0, 2) * pow(sigma_p_pingx, 4) * sigma_p_chuiz * pow(tan(theta), 2) + pow(sigma_p_pingx * sigma_p_chuiz, 2) * sigma_p_chuiz) / pow(pow(sigma_p_pingx * tan(theta), 2) + pow(sigma_p_chuiz, 2), 2.5);
    return exp(-(pow(p_0 * tan(theta), 2) / 2 / (pow(sigma_p_pingx * tan(theta), 2) + pow(sigma_p_chuiz, 2)))) * tan(theta) / pow(cos(theta), 2) / pow(sigma_p_chuiz, 2) / sigma_p_pingx / sqrt(2 * M_PI) * f_theta_integra(0, 100 * sigma_p_chuiz, theta);
}
double data_calculation::px_py_p_chuiz(double p_chuiz)
{
    return p_chuiz / pow(sigma_p_chuiz, 2) * exp(-pow(p_chuiz / sigma_p_chuiz, 2) / 2);
}
void data_calculation::plot_px_py_p_chuiz()
{
    double sum_1 = 0;
    double sum_2 = 0;
    for (double p_chuiz = begin_points; p_chuiz < final_points; p_chuiz += delta_points)
    {
        double f_1 = f_p_jins(p_chuiz);
        // double f_2 = f_theta_2(p_chuiz);
        sum_1 += delta_points * f_1;
        // sum_2 += delta_points * f_2;
        graph->AddPoint(p_chuiz, f_1);
        cout << "p_chuiz:" << p_chuiz << " " << "f(1):" << f_1 << endl;
    }
    cout << sum_1 << endl;
    // cout << sum_2 << endl;
    c1->Divide(1, 1);
    c1->cd(1);
    graph->GetXaxis()->SetTitle("#theta");
    graph->GetYaxis()->SetTitle("f(#theta)");
    graph->Draw("ACP");
    // graph_2->Draw("Same");
}
#endif // !